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We present exact results on the dynamics of a biased, by an external force F, intruder (BI) in 
a two-dimensional lattice gas of unbiased, randomly moving hard-core particles. Going beyond 
the usual analysis of the force-velocity relation, we study the probability distribution P(R n ) of 
the BI displacement R„ at time n. We show that despite the fact that the BI drives the gas to 
a non-equilibrium steady-state, P(R„) converges to a Gaussian distribution as n — > oo. We find 
that the variance a 2 of P(R n ) along F exhibits a weakly superdiffusive growth a 2 ~ i/in ln(n), 
and a usual diffusive growth, a 2 . ~ V2 n, in the perpendicular direction. We determine v\ and 
exactly for arbitrary bias, in the lowest order in the density of vacancies, and show that v\ ~ |F| 2 
for small bias, which signifies that superdiffusive behaviour emerges beyond the linear-response 
approximation. Monte Carlo simulations confirm our analytical results, and reveal a striking field- 
induced superdiffusive behavior a 2 ~ n 3//2 for infinitely long 2D stripes and 3D capillaries. 

PACS numbers: 02.50.-r, 05.40.-a, 66.30.-h, 68.35.Fx 



A biased intruder (BI) traveling through a quies- 
cent medium, composed of bath particles which move 
randomly without any preferential direction, drives the 
medium to a non-equilibrium steady-state with an inho- 
mogencous particles' spatial distribution: the latter jam 
in front of and are depleted behind the BI. The BI can 
be a charge carrier biased by an electric field or a colloid 
moved with an optical tweezer. The bath particles may 
be colloids in a solvent or adatoms on a solid surface. 

Such microstructural changes of the medium (MCM) 
were experimentally observed, e.g., in microrheological 
studies of the drag force on a colloid driven through a A- 
DNA solution [l[ or for a BI in a monolayer of vibrated 
grains Q ■ Brownian Dynamics simulations have revealed 
the MCM for a driven colloid in a A-DNA solution \Mk, 
or for Bis in colloidal crystals Q • Remarkably, the MCM 
not only enhance the drag force exerted on the BI, but 
also induce effective interactions between the Bis, when 
more than one BI is present 

The MCM produced by a BI were studied analytically 
for quiescent baths modeled as hard-core lattice gases 
with symmetric simple exclusion dynamics 0-[l5|. De- 
spite some simplifications (interactions are a mere hard- 
core, no momentum transfer and etc), this type of mod- 
eling captures quite well many qualitative features and 
reproduces a cooperative, essentially many-particle be- 
havior present in realistic physical systems [la . Il7 |. It is 
also often amenable to analytical analysis. 

It was found that in ID the size of the inhomogene- 
ity grows in proportion to the traveled distance, so that 
the jamming- induced contribution to the frictional drag 
force exhibits an unbounded growth with time n. In con- 
sequence, the BI velocity v n oc n -1 / 2 [§-ll|, ensuring 



the validity of the Einstein relation for anomalous tracer 
diffusion in ID hard-core lattice gases 
higher dimensions, the BI velocity approaches a termi- 
nal value v and the bath particles' distribution attains a 
non-equilibrium stationary form with a jammed region in 
front of the BI and a depicted region in its wake [13l4l5| | . 
Strikingly, behind the BI the bath particles density ap- 
proaches the mean value p as a power-law of the distance 
x: 1/rr 3 / 2 and h\(x)/x 2 in 2D and 3D [IMBl, which sig- 
nifies that the medium remembers the passage of the BI 
on large temporal and spatial scales. 

Given such an anisotropy in the bath particles' distri- 
bution, one might be curious if the distribution P(R„) of 
the BI vector displacement R„ — (X, Y) would also be 
asymmetric along the axis of the applied force. In this 
paper we address the question of the asymptotic forms 
of the propagator and provide an exact solution within 
the lattice gas model with a high bath particles density 
p, (or small density po = 1 — p of vacancies), i.e. for an 
asymmetric simple exclusion process (ASEP) evolving in 
a 2D dense sea of symmetric exclusion processes (SEPs). 

Using the anal ytic al approach developed previously by 
two of us in Ref.[l9(, we set out to show that, curiously 
enough, in the lowest order in the density of vacancies 
P(R„) along the axis of the applied force tends towards 
a Gaussian function as n — > oo. Clearly, when there is 
no external bias so that the bath is homogeneous, con- 
vergence to a Gaussian is quite trivial (2fj| . In our case, 
however, given a non-equilibrium situation and essential 
MCM, this result certainly cannot be expected a pri- 
ori. More striking, we will show that the variance ct 2 
of P(R„) in the direction of the applied bias grows at a 
faster rate due to an additional logarithmic factor, than 
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in the direction perpendicular to the force, which is a 
manifestation of a rather counter-intuitive cooperative 
behavior. Monte Carlo simulations confirm our analyt- 
ical results, and reveal a surprising superdiffusive field- 
induced growth cr 2 ~ n 3 / 2 for infinite 2D stripes and 3D 
capillaries. On contrary, our simulations show that for 
single- file geometries cr 2 is completely independent of the 
bias. To the best of our knowledge, these results are new. 

Consider a square lattice of L x x L y sites r = (x, y) 
with integer valued components and periodic boundary 
conditions. The lattice is populated with hard-core bath 
particles and a single hard-core intruder is initially at the 
origin. M lattice sites are vacant. The system evolves in 
discrete time n and particles move randomly (subject to 
hard-core exclusion) by exchanging their positions with 
the vacancies. The bath particles have symmetric hop- 
ping probabilities, i.e., given a vacancy is at an adjacent 
site, any bath particles exchanges its position with the 
vacancy with probability = 1/4 independently of the di- 
rection. On the other hand, the intruder is subject to 
a constant force F oriented in the positive x direction. 
The normalized jump probabilities of the {isolated) BI 
are given, in the usual fashion (see, e.g., Ref.[lj|), by 



p v = exp(^(F-e u )\ /£exp(J(F.e„)) , (1 



where f3 is the reciprocal temperature, e„ is the unit vec- 
tor denoting the jump direction, v € {±x, iy}, (F-e„) is 
a scalar product and F = F e x . The sum with the sub- 
script /i (the normalization constant) denotes summation 
over all possible orientations of the vector e M . 

We turn now to the limit of small density of vacancies, 
po = M/(L X x L y ) <C 1, and focus on the behavior in 
the lowest order in pq. Then, it is expedient to formulate 
the dynamics of the system in terms of the dynamics of 
vacancies. Following Ref.[2(J, we stipulate that at each 
tick of the clock each vacancy makes a step exchang- 
ing its position with a bath particle chosen at random 
(with probability = 1/4) from among its four neighbors, 
in case when neither of them is the BI. In case when one 
of the neighboring particles is the BI, the situation is a 
bit more complicated. According to Ref. [l!| a correct 
choice of the transition probabilities, which avoids spu- 
rious temporal trapping, is as follows: if a vacancy is 
at site R„ + e„ at time moment n and the BI occupies 
site R„, then it exchanges its position with the BI with 
probability c/_„ = p v j (3/4+p v ), and with the probability 
= 1/(3 + 4p„) with any of three adjacent bath particles. 
Note that in a complete description of the lattice gas dy- 
namics, these rules would have to be supplemented for 
cases where two vacancies are adjacent or have common 
neighbors; however, these cases contribute only to O{p v ) 
so that we can leave the rules for such events unstated. 

The steps involved in the derivation of P(R n ) arc 
discussed in Ref. fljj . In the thermodynamic limit, 



(M,L x ,L y — > oo with po <C 1 kept fixed) the normal- 
ized P(R n ) is given by 



P(R n ) ~ J dk exp (-i ( k ' R ») - P0fin(k)) , 

with tt n (k) defined via its Z-transform : 

*(k) 



(i-z) 2 (i + x - 1 *(k)) 



(2) 
(3) 



where " ~ " denotes the exact leading at 2— >l _ (n^>l) 
behavior, Xz ~ — 7r/ln(l — z) is the generating function 
of the mean number \n of " new" sites visited at the n-th 
step (see, e.g., Ref. [5o|]), *(k) = -ia k x + aifc 2 /2 + 
a2fc?,/2 and a,j are given, for arbitrary £ = /3P/2, by 



a 



sinh(£) 



(2n - 3) cosh(£) + 



"2 



- , at = coth(£) a , (4) 



(5) 



cosh(£) + 2?r - 3 ' 



Differentiating $(k) = exp (— po&l n (k)) with respect to 
k x and k y , we find that the BI mean velocity v and the 



variances cr 2 and o~ x along the y- and x-axes obey 



Po a , 
Po 0-2 n , 

f 2 «o ,„ 
Po ai H [Hn+x 

\ 7T 
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where H n = Ylk—i fc 1 is the n-th harmonic number. 

In Fig. [T] (panels a and b) we present a comparison 
of our analytical results for the mean velocity and the 
variance in the y-direction and the results of numerical 
simulations for fixed bias f3F = 5 and different values of 
Pq. One notices that our analytical results are in a very 
good agreement with the numerical data for po up to 0.2. 

Noticing next that in the large n limit H n+ \ ~ ln(n) + 
7 + 0(l/n), where 7 sa 0.577 is the Euler-Maschcroni 
constant, we find that asymptotically 



Po cti + 



M (7 

7T 
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V ln(n) ) 



(9) 



On comparing the results in Eqs. Q and ©, we thus 
conclude that a x grows faster due to an additional log- 
arithmic factor than cr 2 ,, which shows the usual diffusive 
behavior. It means that, rather counter-intuitively, the 
distribution becomes progressively more broad along the 
x direction, compared to its behavior along the y direc- 
tion. Note that, since cio ~ /3F for /3F <C 1, the coef- 
ficient before the superdiffusive term is proportional to 
(f3F) 2 , which signifies that the superdiffusive behavior 
emerges beyond (and thus can not be obtained within) 
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FIG. 1. (color online) Velocity v (panel a) and <Jy/n (panel 6) 
versus po for j3F = 5. The dashed lines define our theoretical 
predictions in Eqs. (|6]) and (0. The symbols are the MC 
simulations results. The inset in panel a is a zoom of the 
region po < 0.01. In panel c we plot (f>(n) (see the text) vs 
n for po = 0.002, fiF = 100 (□), /3F = 5 (<), /3F = 2 (O) 
and pF = 1 (O). The dashed curves are our predictions in 
Eq. @. 



the linear-response approximation. 

In Fig. [1] (panel c) we present an evidence for the addi- 
tional logarithmic factor in the x-component of the vari- 
ance. To single out this contribution, we plot here versus 
ln(n) the function cj)(n) = a x /n — po( a i + 2ap(7 — l)/7r), 
which according to Eq. ((9]) should grow in proportion to 
ln(n). One observes an apparent crossover to a logarith- 
mic behavior which persists then over two time decades. 

The skewncss 71 (x) of P(X) = J dY P(R n ) obeys 



6a (a 2 , + TOi/ hi(rt)) j ln(n) 
(2a 2 + 7tai/ ln(n)) 



3/2 



7rp ™ 



(10) 



Note that 71 (x) > so that P{X) has a positive skew 
which implies that fluctuations in the BI position are 
more pronounced for X > v n (i.e. in the region where 
the bath particles jam) than for X < v n where the bath 
particles are depleted. Next, for the kurtosis we get 



6(4a^ + 6a^7rai/ln(n) +7r 2 af/ln 2 (n)) ln(n) 
(2a 2 , + nai/ ln(n)) 2 



Trpon 
(11) 

i.e., 72(2;) vanishes as n — > 00, and hence, P(X) converges 
to a Gaussian, despite a non-equilibrium and MCM! 

We find next that the skewness 71 (y) of PiY) 
j dX P(R n ) vanishes. For the kurtosis of P(Y) we ob- 
tain 72(2/) <~ 61n(n)/7rpo ri , which implies that 72(2/) de- 
cays exactly in the same way as 72(2:) and moreover, 
J2{y)/j2(x) — > 1~ as n — > 00. Observe that 72(2/) is 
independent of the bias, while 71 (x) and 72(2;) become 
independent of it when ln(n) 3> nai. 

In Fig. ^ (panels a and 6) we compare our analyti- 
cal predictions for 71 (x) and 72 (x) in Eqs. (fTU)) and (fTTj) 
against the results of numerical simulations (symbols). 
To single out the logarithmic factors, we plot here the 
functions 71 (x) = 117^(3;) and 72 (x) = 7172 (x) vs ln(n). 
Note that the kurtosis approaches the asymptotic predic- 
tion in Eq. (|lip quite rapidly, at times of order n ~ 10 2 
(and moreover, the /J.F'-independent behavior is estab- 
lished quite fast as well). This confirms our result that 
the kurtosis vanishes as 72 (x) oc ln(n)/n and, hence, that 
the distribution in the direction of the bias converges to a 
Gaussian. The skewness, which decays at a slower rate, 
approaches the asymptotic result in Eq. (JTUJ) a decade 
later and the /3.F-independcnt behavior sets in at larger 
times, which are not accessible in our numerical simula- 
tions of this essentially many-particle system. 

We turn next to the large- (but finite) n correc- 
tions to the asymptotic Gaussian distribution. Multi- 
plying both denominator and the nominator of Eq. Q 
by the complex conjugate of the denominator, introduc- 
ing a small parameter e x = a\ ln(n)/27rer 2 oc 1/n (e y = 
02 ln(n)/27rcr 2 oc ln(n)/n) and expanding <&(k Xl k y = 0) 
(<&(k x = 0, k y )) up to the second order in e x (e y ), we get 



P(X) = 



exp f — (n x — r) x ) 2 /2 ) . f . . 
V L { 1 + [3 (g - f) + 6 ft - (g + /) t x + 2 (3 (g - 2f) + (g + 2/) rf x ) rj x n x 



(g-f + ft) Vl - 2(5 - 2/) rj x 4 + (g - /) nt] e x + 0{e 2 x ln(n))} 



(12) 



^)- CXP( r ^- /2) {l+[3-6,g + <] ^+0(el)\ 



(13) 
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where we have conveniently chosen n x = X/o~ x {rj y = 
Y / a y)> V x = vn / a x, .9 = 1 + Ooln(n)/7rai and / = 
7rai/2(7rai + 2agln(n)). Note that the result in Eq. (fP2")l 
becomes identical to that in Eq. ([13)) when /3.F = 0. 




FIG. 2. (color online) Reduced skewness 71 (a; ) = nj 2 (x) 
(panel a) and the kurtosis 72(2;) = 7172(2;) (panel b) for po — 
0.002, /3f = 1 (O) and jSF = 100 (□). Dashed and solid 
curves are the corresponding theoretical results in Eqs. (|10[) 
and Panels c and d present the time evolution of the 

integrated distributions: P(X) for p = 0.002 and fiF = 100 
for times n = 10 4 ( ), 2 x 10 4 (v), 3 x 10 4 (+), 4 x 10 4 (O) 
and 10 5 (O). The inset shows the deviation AP(X) of our 
result in Eq. (|12[1 from the numerical data. Panel 6: P{Y) 
for p = 0.002 and /3F = 1 for n = 10 4 ( ), 1.5 x 10 4 (A) 
and 2 x 10 4 (v)- Solid curves are our results in Eqs. (|12p and 
(|13|l . The black dashed line defines the asymptotic Gaussian 
distribution for n — 2 x 10 4 . 

In Fig. [5] (panels c and d) we compare our Eqs. (fT2")) 
and (fT3")) against the numerical simulations data. We ob- 
serve that as time progresses, the discrepancy between 
Eq. (|T2|) and numerically obtained P(X) gets smaller, 
as evidenced by the inset in panel c. The convergence 
is non-uniform so that we have a better agreement be- 
tween the theory and numerics for the values of X to 
the left from the maximum, than to the right from it 
(recall that 7i(x) > 0). Along the y-axis (panel gQ, we 
observe a pretty good agreement between our Eq. (fT3|) 
and the numerical data. Note, however, that some slight 
discrepancy between a Gaussian (blue dashed line) and 
a Gaussian with the correction term, Eq. (|13[) . persists 
even for the longest observation time, n = 2 x 10 4 . 

The following remark is in order: An appearance of 
the logarithmic factor can be interpreted as a sign that 
d = 2 is the marginal dimension for this system, above 
which (e.g., in 3D) one should expect a usual diffusive 



growth of a 2 (with, however, a larger prefactor than in 
the perpendicular to the bias directions), while in ID, i.e., 
for a single-file diffusion, the effect should be stronger 
and a 2 would get an additional power-law dependence 
on time. For 3D, indeed, such a scenario seems plausi- 
ble. We have numerically simulated a BI dynamics in a 
single-file symmetric lattice gas and realized that, sur- 
prisingly, the growth law of a 2 is not affected by the 
bias, i.e., a 2 ~ n 1 / 2 . Moreover, our simulations show 
that even the prefactor in this dependence is independent 
of /3F. On contrary, our preliminary estimates suggest 
that a strongly superdiffusive behaviour a 2 ~ n 3 / 2 will 
take place in quasi-lD systems - infinite two-dimensional 
stripes (L x = 00 and finite L y ) and three-dimensional 
rectangular capillaries (L x = 00 and finite L y and L z ). 
Indeed, our Eq. (j9|) tells that a 2 ~ n/Xn- By defini- 
tion, Xn — S n — S n —i, where S n is the mean number 
of distinct sites visited by any individual vacancy up 
to time n. For infinite stripes one has S n ~ n 1 / 2 /L y , 
while S„ ~ n 1 / 2 / L y L z for infinite rectangular capillar- 
ies. Therefore, in both cases we have a 2 ~ n 3 ^ 2 which 
strongly superdiffusive behaviour is unambiguously con- 
firmed by our numerical simulations (see, Fig. 13]). 

We note finally that a similar superdiffusive behavior 
a 2 ~ rf with Q w 1.45 has been observed in Molecular 
Dynamics simulations of a biased intruder dynamics in 
a glass-for ming binary Yukawa liquid in [2l| and more 
recently, in 1221 Taken together, our results and the re- 
sults of Refs.|2ll. [22l | hint at a possibility of encountering 
a new physical phenomenon - a field-induced superdiffu- 
sive broadening of fluctuations in dynamics of a biased 
intruder in quiescent molecular crowding environments. 
This is rather surprising, since such dense environments 
are commonly associated with a subdiffusive behavior. 
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FIG. 3. (color online) Monte Carlo simulations results for the 
variance of the BI displacement in stripes with L x = 10 4 
and L y = 3 (panel a) and rectangular capillaries with L x = 
10 4 and L y = L 2 = 3 (panel b) for p = 0.002 and /3F = 100. 
The dashed curves in both panels have the slope n 3 ^ 2 . 
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